% \documentclass[12pt]{article}
% \documentstyle{l-aa}
% \input{l-aa.sty}
% \documentclass{aa}
\documentclass[twocolumn]{aa}
% \input epsf
  \usepackage{graphicx}

% from http://stud3.tuwien.ac.at/~e9825610/LaTeXmath.html

\DeclareMathAlphabet{\mathitbf}{OML}{cmm}{b}{it}

\DeclareMathSymbol{\alpha}{\mathalpha}{letters}{"0B}
\DeclareMathSymbol{\beta}{\mathalpha}{letters}{"0C}
\DeclareMathSymbol{\gamma}{\mathalpha}{letters}{"0D}
\DeclareMathSymbol{\delta}{\mathalpha}{letters}{"0E}
\DeclareMathSymbol{\epsilon}{\mathalpha}{letters}{"0F}
\DeclareMathSymbol{\zeta}{\mathalpha}{letters}{"10}
\DeclareMathSymbol{\eta}{\mathalpha}{letters}{"11}
\DeclareMathSymbol{\theta}{\mathalpha}{letters}{"12}
\DeclareMathSymbol{\iota}{\mathalpha}{letters}{"13}
\DeclareMathSymbol{\kappa}{\mathalpha}{letters}{"14}
\DeclareMathSymbol{\lambda}{\mathalpha}{letters}{"15}
\DeclareMathSymbol{\mu}{\mathalpha}{letters}{"16}
\DeclareMathSymbol{\nu}{\mathalpha}{letters}{"17}
\DeclareMathSymbol{\xi}{\mathalpha}{letters}{"18}
\DeclareMathSymbol{\pi}{\mathalpha}{letters}{"19}
\DeclareMathSymbol{\rho}{\mathalpha}{letters}{"1A}
\DeclareMathSymbol{\sigma}{\mathalpha}{letters}{"1B}
\DeclareMathSymbol{\tau}{\mathalpha}{letters}{"1C}
\DeclareMathSymbol{\upsilon}{\mathalpha}{letters}{"1D}
\DeclareMathSymbol{\phi}{\mathalpha}{letters}{"1E}
\DeclareMathSymbol{\chi}{\mathalpha}{letters}{"1F}
\DeclareMathSymbol{\psi}{\mathalpha}{letters}{"20}
\DeclareMathSymbol{\omega}{\mathalpha}{letters}{"21}
\DeclareMathSymbol{\varepsilon}{\mathalpha}{letters}{"22}
\DeclareMathSymbol{\vartheta}{\mathalpha}{letters}{"23}
\DeclareMathSymbol{\varpi}{\mathalpha}{letters}{"24}
\DeclareMathSymbol{\varrho}{\mathalpha}{letters}{"25}
\DeclareMathSymbol{\varsigma}{\mathalpha}{letters}{"26}
\DeclareMathSymbol{\varphi}{\mathalpha}{letters}{"27}

% Since the roman alphabet does not contain roman lowercase greek, applying \mathrm to any of these redefined commands yields unexpected results. The same holds for \mathsf, \mathtt, … (Basically this is the same if you did \mathcal{123})

% You may want to have a look at the package fixmath, which incorporates these changes.
%
% There are some other glyphs that may be useful, too:

\DeclareMathSymbol{\imath}{\mathalpha}{letters}{"7B}
\DeclareMathSymbol{\jmath}{\mathalpha}{letters}{"7C}


\begin{document}
% \thesaurus{
%             (03.13.4; %Methods: numerical
%              08.18.1  %Stars: rotation
%             )
%           }
\title{ A Newton iteration method for obtaining equilibria
        of rapidly rotating stars}
\author{A. G. Aksenov\inst{1,2} \and S. I. Blinnikov\inst{1,3}}
\offprints{A. G. Aksenov}
\institute{
    Institute for Theoretical and Experimental Physics,
    117218, Moscow, Russia
\and
    Max Planck Institut f\"ur Astrophysik,
    D-85740 Garching bei M\"unchen, Germany
\and
    University of California Observatories/Lick Observatory
    Board of Studies in Astronomy and Astrophysics
    University of California, Santa Cruz, CA 95064
}
   \date{Received 23 December 1993 / Accepted 9 May 1994}


%    \begin{abstract}
\abstract
{}
{A new computational method for the construction of
rapidly rotating stellar models is described; it combines
the high (quadratic) convergence
of the Newton-Raphson procedure with the use of a fine mesh,
typical of self-consistent field schemes.}
{For the presentation of the gravitational potential a finite difference
approximation to the Poisson equation is used, and
the iterations of the potential are used to find
the distribution of matter consistent with its gravity.}
{It is necessary to solve a huge system of linear
equations on every iteration step, but the matrix of the system
is sparse, so that powerful sparse matrix solvers can be applied here.
The high convergence rate of the iteration procedure and
a robust sparse system solver warrant the
high accuracy of calculations.}
{The method is powerful enough for the description of models
with a large density contrast and high $T/|W|$, where
$T$ and $W$ are total kinetic and gravitational energy respectively.
The efficiency of the method is tested in calculations of
the stationary states for polytropes for $T/|W|$ up to $0.425$.}
   \keywords{numerical methods -- stars:
             interiors -- stars: rotation
            }
%    \end{abstract}

   \maketitle


\section{Introduction}
%______________________________________________________________________________
Many numerical methods have been proposed for calculation of two- or
three-dimensional self-gravitating equilibrium structures since the
pioneer work done by James (1964).
One of the most advanced, the Self-Consistent Field (SCF) method developed
by Ostriker \& Mark
(1968), adopted an integral representation of the gravitational potential
and used simple iterations for obtaining the equilibrium.
Many useful modifications of this method have been published, e.g. by
Clement (1974), Blinnikov(1975a,b) and Galkin et al. (1993); and the SCF method
has been applied to many problems (e.g.\ Ostriker \&  Bodenheimer 1968,
Bodenheimer 1971,  Bodenheimer \& Ostriker 1973, Bisnovatyi-Kogan \&
Blinnikov 1974).
An interesting version of this method has been provided by Hachisu
(1986a); it makes it possible to obtain solutions with high
$T/|W|>0.25$ (here $T$ and $W$ are total kinetic and gravitational
energy of a star, respectively) or in the case of high density contrast.
By way of this method, referred to as HSCF, some two- and three-dimensional
single-
and multi-body models have been calculated (Hachisu 1986a,b;
Hachisu et al. 1986a,b;
Eriguchi \& M\"uller 1991, Eriguchi \& M\"uller 1993).
Nevertheless, both SCF and HSCF  have slow (linear)
convergence levels in general.

The fast (quadratic) convergence
of the Newton-Raphson method constitutes another class of self-consistent
techniques.
In this class, Stoeckly's (1965) was the first to develop a method
which used ``accurate'' formulae and a few
(i.e., only $6$ in the original work)
points for representing the difference equations in
the direction of the polar angle, $\theta$,
and many points and simple formulae along the radius $r$.
The solution was obtained using the Newton-Raphson iteration procedure
for the potential.
On the one hand, the use of a rough grid in $\theta$-direction significantly
reduces the rank of the linear system obtained, but on the other hand,
the method fails in the case of highly deformed nonspherical structures.

Eriguchi \& M\"uller (1985a) have proposed another universal and
powerful method using the  Newton-Raphson iteration procedure, which
allows the calculation of highly deformed configurations.
For the iterations, density values in grid points are used.
It requires the solution of a  system of linear equations with a very huge
matrix whose rank is equal to the number of grid points but, contrary to
the method  discussed in the present paper,
is dense and not sparse.
The really adopted two-dimensional grid was only $40\times 15$
even when using a supercomputer ($44\times 16$ or $64\times 12$ in other
applications, see Eriguchi \& M\"uller 1985b, and  M\"uller \& Eriguchi 1985).
Hence  even in the case of
two-dimensional problems, this technique does not allow us to obtain
a high--resolution solution,
although the convergence of the method is good.

Recently, considerable progress has been attained in the numerical solution of
linear systems of equations with sparse matrices of coefficients for
a huge rank (some tens of thousands).
This makes it possible to revert to the Newton-Raphson approach to
the problem under discussion.
If the Poisson equation is approximated by finite differences on the grid
(instead of using the integral representation of potential)
and an appropriate iteration procedure is used,
then it is possible to reduce the problem to that of
solving a set of linear equations with a sparse matrix.
In this paper we suggest such a  method.

Our approach is very similar to the Newton method
developed by Clement (1974),
who used a finite-difference approximation for the Poisson equation and
considered both SCF and Newton's iterations, although
his technique could be applied only to the situations where the sparse
matrix had
a special structure. The number of mesh points for the $\theta$
coordinate was small, and the method failed for $T/|W| \ge 0.25$.
 Our method allows the use of the matrix with any
pattern (only its sparsity is important with respect to efficiency), so
that it
can be generalized to solve various physical problems. At the same time
it is also close to HSCF method
but  has a higher convergence rate and it is not restricted by the
number of spherical harmonics in the expansion of the potential.
It allows to obtain axisymmetric solutions with a high level of
accuracy, even on a scalar computer.

Our work is a part of the
project exploring the scenario for the supernova explosion mechanism
suggested by Imshennik (1992). It requires us to carry out
2- and 3-dimensional hydrodynamic calculations and to test the
stability of very rapidly rotating models in this scenario (Aksenov
\& Imshennik 1993, Imshennik \& Blinnikov 1993).
Our technique allows us to construct very rapidly rotating configurations
almost up to the limit $T/|W|=0.5$; moreover, it is generalizable for
3-dimensional problems.
The limiting configurations are needed also to obtain the integration
constant for an accretion disk (Bisnovatyi-Kogan 1993).

\section{Basic equations}
%______________________________________________________________________________
For quantities in regular units we will use math bold italics.
All dimensionless quantities will be in standard math italics and greek.

We consider stationary equilibria of axially-symmetric
rotating self-gravitational gases:
\begin{eqnarray}
   \mathitbf{  \rho (\vec{v}\nabla )\vec{v}
    +\nabla P+\rho\nabla\Phi}=0, \label{eq:stat0}
%      \rho (\mbox{\boldmath$v$}\nabla )\mbox{\boldmath$v$}
%     +\nabla P+\rho\nabla\Phi=0, \label{eq:stat0}
\end{eqnarray}
where $\mathitbf{\rho}$ is the density, $\mathitbf{P}$ is the pressure, {$\vec{v}$} is the
velocity of matter, and the
gravitational potential $\mathitbf{\Phi}$ satisfies the Poisson equation
\begin{equation}
    \mathitbf{\nabla^2 \Phi = 4\pi G\rho}. \label{eq:Poisson}
\end{equation}
Let $( \mathitbf{r}, \theta, \varphi)$ be standard spherical coordinates
and ($ \mathitbf{\varpi, z}, \varphi)$ be cylindrical coordinates.
We assume the fluid velocity without meridional circulation
\begin{equation}
    \mathitbf{ \vec{v}
    =\varpi\Omega(\varpi,z)\vec{e}_{\varphi}},
\end{equation}
(here  $ \mathitbf{\Omega}$ is the angular velocity)
and a barotropic equation of state
\begin{equation}
    \mathitbf{ P=P(\rho)}. \label{eq:bareos}
\end{equation}
In this case it is easy to prove the Poincar\'e theorem (see
 e.g. Tassoul 1978) on the constancy of the angular velocity on
coaxial cylinders:
\begin{equation}
     \mathitbf{\partial\Omega(\varpi,z)/\partial z} = 0 \; . \label{eq:poinc}
\end{equation}
Then one can obtain  the Bernoulli
integral of equation (\ref{eq:stat0})
\begin{equation}
%     H(\rho)+\Phi+\Psi=C, \label{eq:stat1}
   \mathitbf{ H(\rho)+\Phi+\Psi=C}, \label{eq:stat1}
\end{equation}
where $C$ is a constant, while the enthalpy
\begin{equation}
     \mathitbf{H(\rho)=\int^{P(\rho)}\frac{{\rm d} P}{\rho}}, \label{eq:enth}
\end{equation}
and the centrifugal potential
\begin{equation}
     \mathitbf{\Psi=-\int^\varpi \Omega^2(\varpi)
          \varpi{\rm d}\varpi}.
\end{equation}

\section{Numerical method}
%____________________________________________________________________________
\subsection{Iteration procedure}
%____________________________________________________________________________
The choice of the basis of dimensionless variables is important for the
numerical scheme. To calculate real stellar models, when the total mass
$ \mathitbf{M}$ must remain constant along a rotating sequence, $ \mathitbf{M}$
shall be specified
as one of the basic dimensional units. In this case our method has to use
a relation similar to Eq.(31) in the paper by Clement (1974), which is
easily incorporated in our scheme. For the barotropic sequences a simpler
approach is possible when
we adopt the following basic dimension units:
$ \mathitbf{G}$ -- the gravitational constant,
$ \mathitbf{R_e}$ -- the equatorial radius of the star, and
$ \mathitbf{\rho}_0$ -- the value of density in a fixed point labeled by the subscript 0.
Then the unit for $ \mathitbf{\Phi, \; \Psi, \; H}$ is $ \mathitbf{\Phi_u=GR_e^2\rho}_0$.
One can express the pressure as $ \mathitbf{P=GR_e^2\rho}_0^2 P$,
the square of the angular velocity as $ \mathitbf{\Omega^2=G\rho}_0\Omega^2$,
the star mass as $ \mathitbf{M=R_e^3\rho}_0 M$,
the angular momentum as $ \mathitbf{J=G^{1/2}R_e^5\rho}_0^{3/2} J$,
the gravitational energy as $ \mathitbf{W=GR_e^5\rho}_0^2 W$,
and the kinetic energy as $ \mathitbf{T=GR_e^5\rho}_0^2 T$.
Here we denote as $F$ the dimensionless form of a variable $ \mathitbf{F}$.
We have to solve the following  system of equations
(the operator $\nabla^2$ below is assumed in dimensionless variables)
\begin{equation}
     \mathitbf{\rho}_0 \rho= \mathitbf{H}^{-1}[ \mathitbf{\Phi_u}( C- \Phi- \Psi)] \label{eq:stat_d},
\end{equation}
\begin{equation}
    \nabla^2  \Phi = 4\pi  \rho. \label{eq:Poisson_d}
\end{equation}
In the point 0: $ \mathitbf{\rho=\rho}_0$,
\begin{equation}
     \mathitbf{\Phi_u}( C- \Phi_0- \Psi_0)= \mathitbf{H(\rho}_0) \label{eq:c0},
\end{equation}
and in the two fixed boundary points $A$, $B$ where $ \mathitbf{\rho}=0$ (Hachisu 1986a):
\begin{equation}
     \mathitbf{\Phi_u}( C- \Phi_A- \Psi_A)= \mathitbf{H}(0), \label{eq:c1}
\end{equation}
\begin{equation}
     \mathitbf{\Phi_u}( C- \Phi_B- \Psi_B)= \mathitbf{H}(0). \label{eq:c2}
\end{equation}
Of course, one can choose the definition of enthalpy in (\ref{eq:enth})
in the most natural way putting $ \mathitbf{H}(0)=0$, but we retain $ \mathitbf{H}(0)$ in
the following  for the sake of generality.
The point $A$ will be chosen as a point on the star surface in the
equatorial plane
so that $\rho(r>r_A)=0$.

It is convenient to take the function $ \Psi$ in the form
\begin{equation}
     \Psi=C_\Psi \chi( \varpi), \label{eq:Psi}
\end{equation}
where $C_\Psi$ is an unknown constant,
and $ \chi( \varpi)$
is a prescribed function, for example (Eriguchi \& M\"uller 1985a),
\begin{equation}
     \chi( \varpi)
        =\cases{- \varpi^2/2
                & rigid rotation,
                \cr
                -\frac{1}{2}
                 \ln(d^2+ \varpi^2)
                & $v$-constant rotation,
                \cr
                -1/[2(d^2+ \varpi^2)]
                & $j$-constant rotation;
         } \label{eq:rolaw}
\end{equation}
here $d$ is some constant.

From conditions (\ref{eq:c0}),(\ref{eq:c1}),(\ref{eq:c2}), it is possible to
express the constants $C_\Psi$, $ C$, $ \Phi_0$ through the values
of $ \Phi$ in points 0, $A$, $B$:
\begin{equation}
    C_\Psi=\frac{ \Phi_B- \Phi_A}
           { \chi_A- \chi_B},
           \label{eq:C_Psi}
\end{equation}
\begin{equation}
     C=\frac{  \Phi_A+ \Psi_A
                 -( \Phi_0+ \Psi_0)\mathitbf{H}(0)/\mathitbf{H}(\mathitbf{\rho}_0)}
           {1-\mathitbf{H}(0)/\mathitbf{H}(\mathitbf{\rho}_0)}, \label{eq:C}
\end{equation}
\begin{equation}
     \mathitbf{\Phi_u}=\frac{ \mathitbf{H(\rho}_0)}{ C- \Phi_0- \Psi_0}, \label{eq:Phi_0}
\end{equation}
and rewrite equation (\ref{eq:stat_d})
\begin{equation}
     \mathitbf{\rho}_0 \rho
    = \mathitbf{\rho_H}\left[\frac{ \mathitbf{H(\rho}_0)}{ C- \Phi_0- \Psi_0}
            ( C- \Phi- \Psi)\right],
\label{eq:stat_d1}
\end{equation}
where we introduce $\mathitbf{\rho_H(.)} \equiv {\mathitbf H^{-1}(.)} $,
that is the inverse function of $ \mathitbf{H(\rho})$ which outputs density given $H$,
and $ \Psi \equiv C_\Psi \chi$, $C_\Psi$, $ C$
are functions of the unknown potential $ \Phi$.
So we have to find $ \Phi$ from
(\ref{eq:Poisson_d}),(\ref{eq:stat_d1}).

We shall use the Newton iteration procedure for this,
for the $(k+1)$-th iteration step equations
(\ref{eq:Poisson_d}), (\ref{eq:stat_d1}) give:
\begin{eqnarray}
    \nabla^2 \Phi^{k+1}=
    4\pi\Bigl\{ \rho^k
	  +\frac{{ \mathitbf{(\rho_H)}}^\prime}{ \mathitbf{\rho}_0}
           \frac{ \mathitbf{H(\rho}_0)}{ C- \Phi_0- \Psi_0} \nonumber \\
     	   \times
           \Bigl[
	        ( C^{k+1}- C^k)
	       -( \Phi^{k+1}- \Phi^k)
	       -( \Psi^{k+1}- \Psi^k)               \nonumber \\
               -\Bigl( ( C^{k+1}-  C^k)
                      -( \Phi_0^{k+1}- \Phi_0^k)
                      -( \Psi_0^{k+1}- \Psi_0^k)
                \Bigr)                                    \nonumber \\
                \times
                \frac{  C- \Phi- \Psi}  % check if inverse is correct as in AA paper
                {  C- \Phi_0- \Psi_0}
           \Bigr]
    \Bigr\}, \label{eq:iter}
\end{eqnarray}
where the quantities without superscripts are taken at the $k$-th step of
iteration, and
$(  C^{k+1}-  C^k)$,
$( \Psi^{k+1}- \Psi^k)$,
$( \Psi_0^{k+1}- \Psi_0^k)$
should be expressed through
$( \Phi_0^{k+1}- \Phi_0^k)$,
$( \Phi_A^{k+1}- \Phi_A^k)$,
and $( \Phi_B^{k+1}- \Phi_B^k)$
from (\ref{eq:Psi}),(\ref{eq:C_Psi}),(\ref{eq:C}).
So, actually, the unknown function is $ \Phi$.

In the general barotropic case, the order of calculations should be the
following: to set $ \mathitbf{\rho}_0$,
in using the iteration procedure (\ref{eq:iter}) to find $ \Phi$,
and then to find
$ \mathitbf{R_e}$ from the relation $ \mathitbf{GR_e^2\rho}_0= \mathitbf{\Phi_u}$.
Here $ \mathitbf{\Phi_u}$ is found from (\ref{eq:Phi_0}).

In the polytropic case:
$ \mathitbf{P=K\rho}^{1+1/n}$, \newline
$ \mathitbf{H(\rho)}=(1+n) \mathitbf{K\rho}^{1/n}$, there is a free scaling parameter;
for our choice of units it is $ \mathitbf{\rho}_0$. It is clear that $ \mathitbf{\rho}_0$
disappears in equation (\ref{eq:stat_d1}), so that
we can find the set of solutions for different $ \mathitbf{\rho}_0$
by solving one problem with fixed points 0, $A$, $B$.

\subsection{Computational grid}
%____________________________________________________________________________

There are two general ways to obtain a finite representation of
continuous fields (like $\Phi$ and $\rho$). One way is to use  base
functions and truncated series for the fields expanded in this basis
(like Ostriker \& Mark 1968). Another way is to use a direct finite
difference representation (like Clement 1974). It is common practice to combine
finite differences for $r$ and expansions in Legendre polynomials
for $\theta$. The basis function approach has many advantages (see
Hernquist \& Ostriker 1992 for the discussion of this approach in
the context of stellar dynamics). We have gained some experience using this
approach in the SCF method: the use of fundamental solutions by way of the
Laplace equation in spheroidal coordinates allows us to obtain very flattened
configurations
(Blinnikov 1975b, Bisnovatyi-Kogan \& Blinnikov 1981)
In the present work we want to develop a robust
method that produces an equilibrium solution in the same representation
as used in multi-dimensional hydrodynamic calculations --- that is, in the
finite-difference one.

For the solution of equation (\ref{eq:iter}) we follow some ideas
developed by Clement (1974), i.e. we introduce a uniform grid
for the radial coordinate inside a sphere, encircling the bulk of the star
mass,
and we use an inversion to solve the Laplace equation up to infinite $r$.
We use the grid $(  r_i,\theta_j)$ in the spherical coordinate system.
The $r$-mesh coordinate is taken to be:
\begin{equation}
      r_i=\frac{i}{i_A},\mbox{  }
    \Delta   r_i\equiv   r_i-  r_{i-1}={\rm constant},
    \mbox{ at }0\leq i\leq i_A,
\end{equation}
where $i_A$ is the point on the star surface in the equatorial plane,
and
\begin{eqnarray}
    \frac{1}{  r_{i-1}}-\frac{1}{  r_i}
    \equiv h = {\rm constant},
    \mbox{ at }i_A\leq i\leq i_{\rm max},
\end{eqnarray}
$  r_{i_{\rm max}}=\infty$.
So,
\begin{eqnarray}
    h=\frac{1}{  r_{i_A-1}}-\frac{1}{  r_{i_A}}
    =\frac{1}{i_A-1}=\frac{1}{  r_{i_{\rm max}-1}} \nonumber \\
    =\frac{1}{i_{\rm max}-i_A} \; ,
\end{eqnarray}
 and
\begin{equation}
    i_A=(i_{\rm max}+1)/2
\end{equation}
should be taken.

The $\theta$-mesh coordinate can be chosen in the form
(the plane $\theta=\pi/2$ is the plane of symmetry):
\begin{equation}
    \theta_j=\frac{\pi}{2}\frac{j}{j_{\rm max}},
              \mbox{ }0\leq j\leq j_{\rm max}.
\end{equation}

The radial part of operator $\nabla^2$:
\begin{equation}
\nabla^2_r \Phi
 \equiv\frac{1}{  r^2}\frac{\partial}{\partial   r}
       (  r^2\frac{\partial \Phi}{\partial   r})
 =\frac{1}{  r^4}\frac{\partial^2 \Phi}{\partial {(1/  r)}^2}.
\end{equation}
In the region $  r<1 \; (i\leq i_A-1)$ we use the following
approximation of $\nabla^2_r$:
\begin{eqnarray}
     \nabla^2_r \Phi(  r_i)
    \simeq\frac{1}{(\Delta   r)^2  r_i^2}
     \Bigl[ \frac{  r_{i+1}^2+  r_i^2}{2}( \Phi_{i+1}- \Phi_i)
     \nonumber \\
      -\frac{  r_i^2+  r_{i-1}^2}{2}( \Phi_i- \Phi_{i-1})
     \Bigr],
\end{eqnarray}
in the region $1\leq  r \; (i_A\leq i\leq i_{\rm max})$:
\begin{equation}
     \nabla^2_r \Phi(  r_i)
    \simeq\frac{1}{  r_i^4 h^2}
     \Bigl[ \Phi_{i-1}-2 \Phi_i+ \Phi_{i+1}\Bigr].
\end{equation}
The angular part of $\nabla^2$:
\begin{equation}
     \nabla^2_\theta \Phi
    \equiv\frac{1}{  r^2\sin\theta}\frac{\partial}{\partial\theta}
          \left(\sin\theta\frac{\partial \Phi}{\partial\theta}\right)
\end{equation}
is approximated by
\begin{eqnarray}
     \nabla^2_\theta \Phi(  r_i,\theta_j)
    \simeq\frac{1}{  r_i^2\sin\theta_j(\Delta\theta)^2}
     \Bigl[ \frac{\sin\theta_{j+1}+\sin\theta_j}{2}
       (  \Phi_{i,j+1} \nonumber \\
        - \Phi_{i,j})
      -\frac{\sin\theta_j+\sin\theta_{j-1}}{2}
       ( \Phi_{i,j}- \Phi_{i,j-1})
     \Bigr].
\end{eqnarray}

Now we describe our treatment of the boundary conditions.
For $  r\rightarrow 0$ we have
\begin{equation}
\int_0^\pi\nabla^2 \Phi\sin\theta d\theta\rightarrow
3\int\frac{\partial^2 \Phi}{\partial  r^2}\sin\theta d\theta \; ,
\end{equation}
so
\begin{eqnarray}
     4\pi \rho_{0,k}
    =\int_0^\pi\nabla^2 \Phi\sin\theta d\theta\mid_{  r=0} \nonumber \\
    \simeq\frac{6}{(\Delta  r)^2}
     \Bigl( \sum_k\
         \Bigl\{
             \cos\bigl[\max(0,\Delta\theta(k-1/2))\bigr] \nonumber \\
            -\cos\bigl[\min(\pi/2,\Delta\theta(k+1/2))\bigr]
         \Bigr\}
         ( \Phi_{1,k}- \Phi_{0,j})
     \Bigr),
\end{eqnarray}
where we use the notation $ \Phi_{i,j}\equiv \Phi(  r_i,\theta_j)$.
For $  r\rightarrow\infty$ we have $ \Phi=0$:
\begin{equation}
     \Phi_{i_{\rm max},j}=0.
\end{equation}
At $\theta=0$:
\begin{equation}
     \nabla^2_\theta \Phi_{i,0}
    =\frac{2}{  r^2}
     \frac{\partial^2 \Phi}{\partial\theta^2}\mid_{i,0}
    \simeq
    \frac{2}{  r_i^2}
    \frac{2 \Phi_{i,1}-2 \Phi_{i,0}}{(\Delta\theta)^2},
\end{equation}
At $\theta=\pi/2$:
\begin{equation}
     \nabla^2_\theta \Phi_{i,j_{\rm max}}
    =\frac{1}{  r^2}
     \frac{\partial^2 \Phi}{\partial\theta^2}\mid_{i,j_{\rm max}}
    \simeq
    \frac{1}{  r_i^2}
    \frac{2 \Phi_{i,j_{\rm max}-1}
             - 2 \Phi_{i,j_{\rm max}}}{(\Delta\theta)^2}.
\end{equation}

From (\ref{eq:iter}) we have the system of linear equations for
finding $ \Phi_{i,j}$.
If we denote $u_k= \Phi_{i,j}$ where $k=k(i,j)=1+i(j_{\rm max}+1)+j$,
then we have the sparse system
\begin{equation}
    A\mbox{\boldmath$u$}=\mbox{\boldmath$f$} \label{eq:lsys}
\end{equation}
with nonzero elements:

\noindent
$a_{i,i}$,
$a_{i-1,i}$,
$a_{i+1,i}$,
$a_{i-(j_{\rm max}+1),i}$,
$a_{i+(j_{\rm max}+1),i}$,
$a_{i,k_{i_0,j_0}}$,
$a_{i,k_{0,0}}$,
$a_{i,k_{i_A,j_A}}$,
$a_{i,k_{i_B,j_B}}$,
$a_{
     k_{0,j_{\rm max}},
     k_{1,0}\leq\alpha\leq k_{1,j_{\rm max}-1}
   }$,

\noindent
where $1\leq i\leq k_{\rm max}=k_{i_{\rm max},j_{\rm max}}
  = (i_{\rm max}+1)(j_{\rm max}+1)$.

For the solution of the sparse system (\ref{eq:lsys}),
we use the elimination technique developed by Zlatev (see
\O sterby \& Zlatev 1983). We use our code designed according to
their specifications: the choice of the pivots in $LU$ decomposition
follows the Zlatev strategy, and the iterative refinement may be used for the
improvement of the solution. We find that the iterative refinement is not
needed
for the value of the Zlatev threshold $\le 10^{-8}$. For
the threshold $\approx 10^{-4}$, the iterative refinement is necessary,
but then the $LU$-decomposed matrix is not so dense as for smaller
thresholds, and the algorithm is $\approx$ 30\% faster.

The elements of the matrix $A$ in (\ref{eq:lsys}) change on every
iteration step, but the structure of $A$ (the position of nonzero
elements) remains the same.
This fact is taken into account by the linear equation solver used.

\section{Numerical example: rotating polytropes}
%______________________________________________________________________________
% Two column table
%------------------------------------------------------------------------------
\begin{table*}
    \caption{Polytropes n=1.5 (rigid rotation, spheroidal configurations)}
    \label{Tab_Ririd}
%    \picplace{3.5cm}
\begin{center}
\begin{tabular}{|c l l l l l l l c|} \hline
$  r_B$ & $C_\Psi$ & $  M$ & $  V$ & $  J$ &
$-  W$ & $  P_{\rm max}$ & $   T/(-  W)$ &  $VT$ \\
\hline
$1.000$ & $0.0000$ & $0.6993$ & $4.16$ & $0.000$ &
$0.4192$ & $0.3765$ & $0.0000$ & $9.4\cdot 10^{-6}$ \\

$0.900$ & $0.1033$ & $0.5974$ & $3.76$ & $0.03762$ &
$0.3207$ & $0.3276$ & $0.01885$ & $7.1\cdot 10^{-6}$ \\

$0.800$ & $0.1913$ & $0.4883$ & $3.28$ & $0.03905$ &
$0.2275$ & $0.2765$ & $0.03754$ & $9.3\cdot 10^{-6}$ \\

$0.750$ & $0.2266$ & $0.4304$ & $3.03$ & $0.03557$ &
$0.1835$ & $0.2499$ & $0.04614$ & $9.6\cdot 10^{-6}$ \\

$0.700$ & $0.2539$ & $0.3700$ & $2.75$ & $0.03013$ &
$0.1421$ & $0.2227$ & $0.05343$ & $9.6\cdot 10^{-6}$ \\

$0.650$ & $0.2705$ & $0.3072$ & $2.41$ & $0.02330$ &
$0.1039$ & $0.1949$ & $0.05832$ & $9.6\cdot 10^{-6}$ \\

$0.625$ & $0.2738$ & $0.2752$ & $2.21$ & $0.01961$ &
$0.08643$ &$0.1807$ & $0.05937$ & $9.1\cdot 10^{-6}$ \\
\hline
\end{tabular}
\end{center}
\end{table*}

In our example, rigidly rotating polytropes with $n=1.5$
with spheroidal structure were calculated.
We adopted $r_0=0$; $r_A=1$, $\theta_A=\pi/2$, and the
boundary point $B$ was taken on the polar axis: $r_B\leq 1$, $\theta_B=0$
(following Hachisu 1986a).

The calculations started from a spherical configuration:
$r_B=1$.
After calculating one model, we reduced $r_B$ and calculated the new model
with $ \Phi$ from the preceding step as an initial  aproximation
for $ \Phi$. The initial approximation for the spherical configuration was
rather arbitrary, for example
\begin{equation}
    \Phi=\cases{
                  -1.5, & $  r<1$,
                  \cr
                  -1,   & $  r=1$,
                  \cr
                  -0.5, & $  r>1.$
                }
\end{equation}
The iterations were continued until the correction for $C$ was small:
$\delta  C < 10^{-7}$.
The number of iterations depended, of course, on the proximity of
the initial approximation to the final solution. In the tables below
we show all actual sequences, so that initial approximations appear rather
crude. Nevertheless, the number of iterations ($\leq 4$) was appreciably lower
than in SCF methods, especially for very high $T/W$.

In Table~\ref{Tab_Ririd} we present the dimensionless values for
$  r_B$,
$C_\Psi$, and for the following integrals:
\begin{equation}
  M=\sum \rho_{ij}\Delta V_{ij}  \; ,
\end{equation}
\begin{equation}
  V=\sum_{ \rho_{ij}>0}\Delta V_{ij}  \; ,
\end{equation}
\begin{equation}
  J=\sum \rho_{ij}  v_{\phi_{ij}}  r_i\sin\theta_j
\Delta V_{ij}  \; ,
\end{equation}
\begin{equation}
-  W=-\frac{1}{2}\sum \rho_{ij} \Phi_{ij}\Delta V_{ij} \; ,
\end{equation}
\begin{equation}
   T
 =\frac{1}{2}\sum \rho_{ij}
                  {v}^2_{\phi_{ij}}\Delta V_{ij}  \; ,
\end{equation}
where
\begin{eqnarray}
 \Delta V_{ij}
 =4\pi
  \frac{ {(  r_i+  r_{\min(i+1,i_{\rm max})})}^3
        -{(  r_{\max(i-1,0)}+  r_i)}^3
  }
  {3\times 2^3}  \nonumber \\
  \times
  \Bigl( \cos\frac{\theta_{\max(j-1,0)}+\theta_j}{2}
   -\cos\frac{\theta_j+\theta_{\min(j+1,j_{\rm max})}}{2}
  \Bigr)  \; ,
\end{eqnarray}
\begin{equation}
  v_\phi= \Omega \varpi
            =\Bigl(-\frac{\partial \Psi}
           {\partial \varpi} \varpi\Bigr)^{1/2}  \; ,
\end{equation}

and for the quantity:
\begin{equation}
  P_{\rm max}=\max_{ij}  P_{ij}  \; ,
\end{equation}

The virial test
\begin{equation}
VT=\frac{|2  T+  W +3 \Pi|}{|  W|} \; ,
\end{equation}
where
\begin{equation}
 \Pi=\sum  P_{ij}\Delta V_{ij}  \; ,
\end{equation}
shows the accuracy of the numerical solution.
On the computational grid used,  $i_{\rm max}=401$,\mbox{ }$j_{\rm max}=50$,
we have $VT<10^{-5}$.
Our results (Table~\ref{Tab_Ririd}) are very close to those obtained
by Hachisu (1986a).

We ran the polytrope $n = 1.5$ for the sequence with differential
rotation (the case of $v$-constant rotation in \ref{eq:rolaw}).
In Fig.~\ref{Fig_v-const} we present  the isopycnic contours for one of
the models
in this sequence with its ``dumbbell'' shape for the axis ratio 0.167.
% One columm figure
%------------------------------------------------------------------------------
 \begin{figure}
   \resizebox{\hsize}{!}{\includegraphics[]{1staf.ps}}
   \caption{The density contours for the $v$-constant law
        for polytrope $n=1.5$ and axis ratio $r_B=0.167$
        ($\log_{10}\rho_{min}=-3$,
         $\log_{10}\rho_{max}=0$,
         $\Delta\log_{10}\rho=0.3$).}
    \label{Fig_v-const}
 \end{figure}
% \begin{figure}
% %    \picplace{6.4cm}
% \epsfxsize=250pt \epsfbox{1staf.ps}
%     \caption{
%         The density contours for the $v$-constant law
%         for polytrope $n=1.5$ and axis ratio $r_B=0.167$
%         ($\log_{10}\rho_{min}=-3$,
%          $\log_{10}\rho_{max}=0$,
%          $\Delta\log_{10}\rho=0.3$).
%     }
%     \label{Fig_v-const}
% \end{figure}

One can compare our results for this model with those of
Hachisu (1986b) in Table 2.
Hachisu used a somewhat different system of units,
based on the maximum value of density $\rho_{\rm max}$, instead of
our $\rho_0$. For the model with $r_B=0.167$, we obtain the following
results in the Hachisu units: $M=0.6129$, $V=2.321$, $J=1.506$,
$T=0.06928$, $W=-0.3341$, $P_{\rm max}=0.1652$. The virial test
for this model is $VT=3.1\cdot10^{-5}$.

The method also allows one to obtain the toroidal structures
(we have done this,
but do not present the results, which are again close to Hachisu's, 1986a).
For this, it is necessary to set the point $B$ on the equatorial plane,
to select the point 0 inside the star (for example, between points $A$, $B$),
and to take an appropriate initial approximation for $ \Phi$.

Using the grid $400\times 51$ the solution requires about 70--100 seconds per
iteration on the Silicon Graphics IRIS-4D workstation and on
the Convex 220 computer (without vectorization). The DEC-Alpha workstation
is 3--4 times faster.

\section{Very rapidly rotating polytropes: self gravitating
   thick and slim disks}
%______________________________________________________________________________
% Two column table
%------------------------------------------------------------------------------
\begin{table*}
    \caption{Rapidly rotating polytropes $n=1.5$
        (spheroidal configurations) for the $\alpha$-constant law
    }
    \label{Tab_Rapid1_5}
%    \picplace{15cm}
\begin{center}
\begin{tabular}{|c l l l l l l l c|} \hline
$  r_B$ & $\alpha$ & $  M$ & $  V$ & $  J$ &
$-  W$ & $  P_{\rm max}$ & $   T/(-  W)$ &  $VT$  \\
\hline
  $1.00$ &   $0.000$ &  $0.700$ &  $4.13$ &  $0.000$ &
  $0.420$ &   $0.377$ &  $0.000$ & $3.8\cdot 10^{-5}$ \\
  $0.80$ &   $0.120$ &  $0.579$ &  $3.54$ &  $0.0581$ &
  $0.304$ &   $0.300$ & $0.0483$ & $2.8\cdot 10^{-5}$ \\
  $0.70$ &   $0.187$ &  $0.516$ &  $3.20$ &  $0.0629$ &
  $0.249$ &   $0.259$ & $0.0772$ & $2.5\cdot 10^{-5}$ \\
  $0.60$ &   $0.262$ &  $0.450$ &  $2.86$ &  $0.0625$ &
  $0.197$ &   $0.217$ &  $0.110$ & $2.0\cdot 10^{-5}$ \\
  $0.50$ &   $0.344$ &  $0.383$ &  $2.48$ &  $0.0581$ &
  $0.148$ &   $0.174$ &  $0.149$ & $1.5\cdot 10^{-5}$ \\
  $0.40$ &   $0.436$ &  $0.315$ &  $2.11$ &  $0.0504$ &
  $0.104$ &   $0.130$ &  $0.194$ & $8.1\cdot 10^{-6}$ \\
  $0.30$ &   $0.542$ &  $0.244$ &  $1.68$ &  $0.0399$ &
 $0.0648$ &  $0.0871$ &  $0.248$ & $5.0\cdot 10^{-6}$ \\
  $0.25$ &   $0.601$ &  $0.207$ &  $1.46$ &  $0.0336$ &
 $0.0479$ &  $0.0665$ &  $0.279$ & $1.5\cdot 10^{-5}$ \\
  $0.20$ &   $0.666$ &  $0.169$ &  $1.22$ &  $0.0269$ &
 $0.0328$ &  $0.0471$ &  $0.314$ & $3.4\cdot 10^{-5}$ \\
  $0.15$ &   $0.737$ &  $0.130$ & $0.963$ &  $0.0196$ &
 $0.0199$ &  $0.0295$ &  $0.353$ & $6.7\cdot 10^{-5}$ \\
  $0.12$ &   $0.783$ &  $0.106$ & $0.800$ &  $0.0151$ &
 $0.0134$ &  $0.0201$ &  $0.379$ & $1.1\cdot 10^{-4}$ \\
  $0.10$ &   $0.816$ & $0.0897$ & $0.689$ &  $0.0121$ &
$0.00963$ &  $0.0146$ &  $0.397$ & $1.6\cdot 10^{-4}$ \\
  $0.09$ &   $0.832$ & $0.0812$ & $0.625$ &  $0.0106$ &
$0.00794$ &  $0.0121$ &  $0.406$ & $2.0\cdot 10^{-4}$ \\
   $0.08$ &   $0.849$ & $0.0726$ & $0.569$ & $0.00911$ &
$0.00638$ & $0.00979$ &  $0.416$ & $2.9\cdot 10^{-4}$ \\
   $0.07$ &   $0.867$ & $0.0639$ & $0.499$ & $0.00764$ &
$0.00496$ & $0.00766$ &  $0.425$ & $4.4\cdot 10^{-4}$ \\
\hline
\end{tabular}
\end{center}
\end{table*}
% Two column table
%------------------------------------------------------------------------------
\begin{table*}
    \caption{Rapidly rotating polytropes $n=3$
        (spheroidal configurations) for the $\alpha$-constant law}
    \label{Tab_Rapid3}
%    \picplace{7cm}
\begin{center}
\begin{tabular}{|c l l l l l l l c|} \hline
$  r_B$ & $\alpha$ & $  M$ & $  V$ & $  J$ &
$-  W$ & $  P_{\rm max}$ & $   T/(-  W)$ &  $VT$  \\
\hline
  $1.00$ &   $0.0000$ &  $0.0775$ &  $4.13$ &  $0.00000$ &
  $0.00901$ &   $0.0662$ &  $0.000$ & $5.0\cdot 10^{-5}$ \\
  $0.80$ &   $0.0639$ &  $0.0661$ &  $3.81$ &  $0.00127$ &
  $0.00686$ &   $0.0562$ & $0.0272$ & $6.0\cdot 10^{-5}$ \\
  $0.70$ &   $0.104$ &  $0.0596$ &  $3.61$ &  $0.00140$ &
  $0.00574$ &   $0.0504$ & $0.0447$ & $6.9\cdot 10^{-5}$ \\
  $0.60$ &   $0.151$ &  $0.0526$ &  $3.38$ &  $0.00142$ &
  $0.00463$ &   $0.0441$ &  $0.0658$ & $8.3\cdot 10^{-5}$ \\
  $0.50$ &   $0.207$ &  $0.0451$ &  $3.13$ &  $0.00134$ &
  $0.00354$ &   $0.0372$ &  $0.0919$ & $1.1\cdot 10^{-4}$ \\
  $0.40$ &   $0.277$ &  $0.0370$ &  $2.84$ &  $0.00118$ &
  $0.00250$ &   $0.0297$ &  $0.125$ & $1.5\cdot 10^{-4}$ \\
  $0.30$ &   $0.366$ &  $0.0286$ &  $2.50$ &  $0.00094$ &
  $0.00158$ &  $0.0216$ &  $0.169$ & $2.2\cdot 10^{-4}$ \\
  $0.25$ &   $0.422$ &  $0.0243$ &  $2.30$ &  $0.00080$ &
  $0.00118$ &  $0.0175$ &  $0.197$ & $2.9\cdot 10^{-4}$ \\
  $0.20$ &   $0.489$ &  $0.0199$ &  $2.08$ &  $0.00066$ &
  $0.00082$ &  $0.0133$ &  $0.231$ & $3.9\cdot 10^{-4}$ \\
  $0.15$ &   $0.570$ &  $0.0156$ & $1.80$ &  $0.00050$ &
  $0.00052$ &  $0.0091$ &  $0.274$ & $5.5\cdot 10^{-4}$ \\
  $0.12$ &   $0.630$ &  $0.0129$ & $1.61$ &  $0.00040$ &
  $0.00036$ &  $0.0067$ &  $0.306$ & $7.0\cdot 10^{-4}$ \\
  $0.10$ &   $0.676$ & $0.0111$ & $1.46$ &  $0.00034$ &
  $0.00027$ &  $0.0051$ &  $0.330$ & $8.4\cdot 10^{-4}$ \\
  $0.09$ &   $0.700$ & $0.0102$ & $1.36$ &  $0.00031$ &
  $0.00023$ &  $0.0044$ &  $0.344$ & $9.2\cdot 10^{-4}$ \\
  $0.08$ &   $0.727$ & $0.0093$ & $1.28$ & $0.00028$ &
  $0.00019$ & $0.0036$ &  $0.358$ & $1.0\cdot 10^{-3}$ \\
  $0.07$ &   $0.755$ & $0.0084$ & $1.19$ & $0.00024$ &
  $0.00016$ & $0.0029$ &  $0.373$ & $1.1\cdot 10^{-3}$ \\
\hline
\end{tabular}
\end{center}
\end{table*}

The described method allows us to construct the rapidly rotation models
very close to the limit
\begin{equation}
    T/|W|=0.5,
\end{equation}
which follows from the virial theorem $2T+W+3\Pi=0$. In order to construct
a sequence of very flattened models, we suggest that the rotation
law be defined in such a form that the centrifugal force forms certain fraction
$\alpha$ of the gravitational force at the equatorial plane
($\theta=\pi/2$):
\begin{equation}
     \alpha\frac{\partial \Phi}{\partial  r}= \Omega^2  r
    =-\frac{\partial \Psi}{\partial  r}, \label{eq:plane}
\end{equation}
where $0\leq\alpha= {\rm constant}<1$.
Let us refer to this as the $\alpha$-constant rotation law.
We can take the centrifugal potential $ \Psi( \varpi)$ from
(\ref{eq:Psi}) in the form:
\begin{equation}
      \Psi( \varpi)
%    =-\int_0^{ \varpi} \Omega^2 rdr
    =-\alpha
     \cases{
          \Phi(  r= \varpi,\pi/2)- \Phi_0,
         & $ \varpi<1$,
     \cr
          \Phi(1,\pi/2)- \Phi_0, & $ \varpi>1$.
     }
\end{equation}
For spheroidal configurations (and for the same fixed points 0, $A$, $B$
as in the preceding section), the conditions
(\ref{eq:c0}),(\ref{eq:c1}),(\ref{eq:c2})
give:
\begin{equation}
      C=\frac{  \Phi_B
                 -\mathitbf{H}(0)/\mathitbf{H}(\mathitbf{\rho}_0) \Phi_0}
           {1-\mathitbf{H}(0)/\mathitbf{H}(\mathitbf{\rho}_0)},
\end{equation}
\begin{equation}
    \mathitbf{\Phi_u}=\frac{\mathitbf{H}(\mathitbf{\rho}_0)}{  C- \Phi_0},
\end{equation}
and
\begin{equation}
      \Psi( \varpi)
     =\cases{
         \frac{ \Phi_B- \Phi_A}{ \Phi_A- \Phi_0}
         ( \Phi(  r= \varpi,\pi/2)- \Phi_0),
         & $ \varpi<1$,
     \cr
          \Phi_B- \Phi_A, & $ \varpi>1$.
     }
\end{equation}

In the main iteration relation (\ref{eq:iter}),
we now have more the complicated expression for the difference
$( \Psi^{k+1}- \Psi^k)$ through the values of $\Phi$.

We find
\begin{eqnarray}
( \Psi^{k+1}- \Psi^k)
= \frac{\partial \Psi}{\partial \Phi_0}
      \left( \Phi_0^{k+1}- \Phi_0^k\right) \nonumber \\
 +\frac{\partial \Psi}{\partial \Phi_A}
      \left( \Phi_A^{k+1}- \Phi_A^k\right)
 +\frac{\partial \Psi}{\partial \Phi_B}
      \left( \Phi_B^{k+1}- \Phi_B^k\right) \nonumber \\
 +\frac{\partial \Psi}
  {\partial \Phi( \varpi,\frac{\pi}{2})}
  \left(  \Phi^{k+1}( \varpi,\frac{\pi}{2})
   - \Phi^k( \varpi,\frac{\pi}{2})\right),
\end{eqnarray}
and $ \Psi_0\equiv 0$.

In the computational grid, we define
$ \varpi_{ij}=  r_i\sin\theta_j$,
and we use the  following approximation:
\begin{equation}
    \Phi(  r={ \varpi}_{ij},\theta=\pi/2)
   \simeq
    c_1 \Phi_{i_{\rm e}-1,j_{\rm max}}
    +c_2 \Phi_{i_{\rm e},j_{\rm max}},
\end{equation}
where $i_{\rm e}$ is the nearest to $ \varpi_{ij}$
point such that
$  r_{i_{\rm e}-1}\leq \varpi_{ij}\leq  r_{i_{\rm e}}$, and
%$c_1=\frac{  r_{i_{\rm e}}- \varpi_{ij}}
 %    {  r_{i_{\rm e}}-  r_{i_{\rm e}-1}}$,
%$c_2=\frac{ \varpi_{ij}-  r_{i_{\rm e}-1}}
 %    {  r_{i_{\rm e}}-  r_{i_{\rm e}-1}}$.
$c_1=(  r_{i_{\rm e}}- \varpi_{ij})/
     (  r_{i_{\rm e}}-  r_{i_{\rm e}-1})$,
$c_2=( \varpi_{ij}-  r_{i_{\rm e}-1})/
     (  r_{i_{\rm e}}-  r_{i_{\rm e}-1})$.
Hence, in the matrix $A$ (\ref{eq:lsys}),
new nonzero elements
$a_{k_{i,j},k_{i_{\rm e}(i,j)-1,j_{\rm max}}}$,
$a_{k_{i,j},k_{i_{\rm e}(i,j),j_{\rm max}}}$,
$0\leq i\leq i_{\rm max},\mbox{ }0\leq j\leq j_{\rm max}$
appear.

In Tables~\ref{Tab_Rapid1_5}, \ref{Tab_Rapid3}
we present our results obtained for the $\alpha$-constant
law for the polytropic indexes $n=1.5$ and $n=3$ on the grid
$200\times 51$.
The notation is the same as for Table~\ref{Tab_Ririd}.
Fig.~\ref{Fig_alpha-const} shows the contours
of isopycnic surfaces for the model with axis ratio $r_B=0.07$ in this
sequence. Unlike in the case of $v$- or $j$-constant laws, the shape of the
configuration for the $\alpha$-constant law is much flatter.
% One columm figure
%------------------------------------------------------------------------------
 \begin{figure}
   \resizebox{\hsize}{!}{\includegraphics[]{2staf.ps}}
   \caption{The density contours for the $\alpha$-constant law
        for polytrope $n=1.5$ and axis ratio $r_B=0.07$
        ($\log_{10}\rho_{min}=-3$,
         $\log_{10}\rho_{max}=0$,
         $\Delta\log_{10}\rho=0.2$).}
    \label{Fig_alpha-const}
 \end{figure}

% \begin{figure}
% %    \picplace{2.7cm}
% \epsfxsize=250pt \epsfbox{2staf.ps}
%     \caption{
%         The density contours for the $\alpha$-constant law
%         for polytrope $n=1.5$ and axis ratio $r_B=0.07$
%         ($\log_{10}\rho_{min}=-3$,
%          $\log_{10}\rho_{max}=0$,
%          $\Delta\log_{10}\rho=0.2$).
%     }
%     \label{Fig_alpha-const}
% \end{figure}

\section{Discussion}
%______________________________________________________________________________

The method has the second order accuracy (relative to the size
of the space zones).
The dependence of solution accuracy ($-  C$, $VT$)
on the number of zones (i.e. on $i_{\rm max}$) is presented in
Table~\ref{accuracy}
for the one-dimensional case (i.e. without rotation) for $n=1.5$.
% One column table
%------------------------------------------------------------------------------
\begin{table}
    \caption{
        Dependence of the algorithm accuracy on the number of grid
        points (one-dimensional case)
    }
    \label{accuracy}
\begin{center}
\begin{tabular}{ | r l l |}
\hline
$i_A$ & $-  C$ & $VT$ \\
\hline
$  10$ & $0.705$ & $5.6\cdot 10^{-3}$ \\
$  20$ & $0.703$ & $1.2\cdot 10^{-3}$ \\
$  50$ & $0.7002$ & $1.6\cdot 10^{-4}$ \\
$ 100$ & $0.69954$ & $3.8\cdot 10^{-5}$ \\
$ 200$ & $0.69932$ & $9.1\cdot 10^{-6}$ \\
$ 500$ & $0.69924$ & $1.4\cdot 10^{-6}$ \\
$1000$ & $0.699221$ & $3.5\cdot 10^{-7}$ \\
$2000$ & $0.699217$ & $8.6\cdot 10^{-8}$ \\
$5000$ & $0.6992153$ & $1.4\cdot 10^{-8}$ \\
\hline
\end{tabular}
\end{center}
\end{table}

One can see that
$VT\propto i_{\rm max}^{-2}\propto(\Delta  r)^2$,
where $\Delta  r$ is the radial zone size.
Thus, for solutions deviating considerably  from the spherical shape,
the grid should satisfy the following condition:
$  r_i\Delta\theta_j\simeq\Delta  r_i$.

The number of iterations is virtually independent of the grid size,
and the CPU time spent on one iteration is  $\propto N_z^\beta$, where
$N_z$ is the number of nonzero elements in the matrix $A$ from (\ref{eq:lsys})
and $\beta<2$.

Though our method is faster than the Newton iteration method developed by
Eriguchi \& M\"uller (1985a), it is still rather slow on large grids.
One way to speed up the algorithm would be to use a vectorized version
of our sparse matrix solver. Our experience shows that, in general,
the use of special built-in routines like {\tt SCATTER} and {\tt GATHER}
allows to save about 50\% of the CPU time on a vector computer.
Another way could be to use other sparse matrix solvers,
designed especially for use on the vector machines. More numerical experiments
are needed here however.

In principle, it is possible to extend this method to the three-dimensional
case (cf. Hachisu 1986b
and Hachisu et  al. 1986 a,b),
but the grid (using our matrix solver and workstations of the current
power) will not be very fine, e.g., $100\times 30\times 30$.

In conclusion we discuss briefly the prospects of application of our
method to construction of the models which are more realistic than
polytropes. The application of our method is straightforward when the
barotropic equation of state (\ref{eq:bareos}) is the good approximation
to reality (e.g. in the case of degenerate stars). For a general equation
of state $P = P(\rho, S)$, where $S$ is the entropy, or  $P = P(\rho, T_p)$,
where $T_p$ is the temperature,
the rotation law (\ref{eq:poinc}) is distinguished, since it
is necessary for the secular stability of a rotating star
(Goldreich \& Schubert 1967, Fricke 1968). Then the surfaces of constant
pressure, density and temperature all coincide and for this pseudobarotropic
situation our method may be generalized following the pioneering work
of Mark (1968),  Jackson (1970)  and Bodenheimer (1971) who considered the
aplication of the SCF method to
real hot stars with radiation. Recently, Eriguchi \& M\"uller (1991, 1993)
have developed the generalization of HSCF method for radiating rapidly
rotating pseudobarotropic stars and toroids with account of meridional circulation.
Such a generalization of our method is also possible. We note a very
important paper by Pavlov \& Yakovlev (1978) on the meridional circulation
for very rapidly rotating stars: their approach together with that of
Eriguchi \& M\"uller (1991) may be combined with our method.
This generalization, as well as application of the method to baroclinic
stars (cf. Ury\=u \& Eriguchi 1994), requires, of course, much additional
work.




\begin{acknowledgements}

This work was stimulated by the discussions with G.Bisnovatyi-Kogan,
V.Im\-shennik, E.M\"uller  and D.Nadyozhin. The code was developed
and the draft version of this paper was written
during the A.A.'s stay at MPA, Garching, and S.B.'s stay at
UCO/Lick observatory, Santa Cruz. We are grateful to
W.Hillebrandt and S.Woosley for their hospitality and for providing excellent
working conditions at their institutions. We thank E.M\"uller
for careful reading and improving of the manuscript, we thank also
M.Murzina and O.Bartunov for their help, and D.Koo for allowing us
to use his DEC-Alpha workstation.
The work is supported by the Bundesforschungsministerium, FRG,
in agreement with the Russian Atomic Energy Ministry, and by the
grants from the California Space Institute (CS-58-93),  the
National Science Foundation (AST-91-15367),
and, in part, from the Russian Foundation for
Fundamental Research  (93-02-17114).
\end{acknowledgements}

\begin{thebibliography}{}

\item[] Aksenov A.G., Imshennik V.S., 1994,
Pis'ma AZh, 20, 32--49 \hfill
(Astronomy Lett.,  20, 24--39)

\item[] Bisnovatyi-Kogan G.S., 1993, A\&A  274, 796-806

\item[] Bisnovatyi-Kogan G.S., Blinnikov S.I., 1974,
    A\&A  31, 391-404

\item[] Bisnovatyi-Kogan G.S., Blinnikov S.I., 1981,
    AZh  58, 312-323  (SvA  25, 175)

\item[] Blinnikov S.I., 1975a, AZh 52, 243-254
    (SvA 19, 151-156)

\item[] Blinnikov S.I., 1975b, Dissertation, Sternberg Astronomical
           Institute, Moscow

   \item[] Bodenheimer P., 1971, ApJ  167, 153--163

   \item[] Bodenheimer P., Ostriker J.P., 1973,
      ApJ 180, 159--169

   \item[] Clement M.J., 1974, ApJ 194, 709--714

   \item[] Eriguchi Y., M\"uller E., 1985a,
      A\&A 146, 260--268

   \item[] Eriguchi Y., M\"uller E., 1985b,
      A\&A 147, 161--168

   \item[] Eriguchi Y., M\"uller E., 1991,
         A\&A  248, 435--447

   \item[] Eriguchi Y., M\"uller E., 1993,
         ApJ  416, 666--678

   \item[] Fricke K.J., 1968, Z.Astroph. 68, 317

   \item[] Galkin S.A., Denissov A.A., Drozdov V.V.,

        Drozdova O.M., 1993, A\&A  269, 255--266

   \item[] Goldreich P., Schubert G., 1967, ApJ 150, 571


   \item[] Hachisu I., 1986a,
      ApJ Supplement Series 61, 479--507

   \item[] Hachisu I., 1986b,
      ApJ Supplement Series 62, 461

   \item[] Hachisu I., Eriguchi Y., Nomoto K., 1986a,
      ApJ 308, 161

   \item[] Hachisu I., Eriguchi Y., Nomoto K., 1986b,
      ApJ 311, 214

   \item[] Hernquist L., Ostriker J.P., 1992, ApJ 386, 375-397

   \item[] Imshennik V.S., 1992, Pis'ma AZh 18, 489--504
         (SvA Lett. 18, 194-201)

   \item[] Imshennik V.S., Blinnikov S.I , 1993, in R.McCray, ed.,
     Proc. of IAU Coll. 145, ``Supernovae'', held in Xian, China

   \item[] Ipser J.R., Managan R.A., 1981,  ApJ 250, 362

   \item[] Jackson S., 1970, ApJ 161, 579-585

   \item[] James R. A., 1964, ApJ 140, 552

   \item[] Mark J.W.-K., 1968, ApJ 154, 627-643

   \item[]  M\"uller E., Eriguchi Y., 1985,
      A\&A 152, 325--335

   \item[]  \O sterby O. \& Zlatev Z., 1983,
      Direct Methods for Sparse Matrices.
      Lecture Notes in Computer Science, Vol. 157,
      Springer, Berlin-Heidelberg-New York-Tokyo

   \item[] Ostriker J.P., Bodenheimer P., 1968,
      ApJ 151, 1089--1098


   \item[] Ostriker J.P., Mark J. W.-K., 1968, ApJ 151, 1075

   \item[] Pavlov G.G., Yakovlev D.G., 1978, AZh 55, 1043--1056
        (SvA 22, 595)

   \item[] Stoeckly R., 1965, ApJ 142, 208

   \item[] Tassoul J.-L., 1978, Theory of Rotating Stars,
      Princeton University Press

   \item[]  Ury\=u K., Eriguchi Y.,  1994, preprint, see MNRAS (1995) 277, 1411

% @ARTICLE{1995MNRAS.277.1411U,
%    author = {{Uryu}, K. and {Eriguchi}, Y.},
%     title = "{Structures of rapidly rotating baroclinic stars - II. an extended numerical method for realistic stellar models}",
%   journal = {\mnras},
%  keywords = {METHODS: NUMERICAL, STARS: ROTATION},
%      year = 1995,
%     month = dec,
%    volume = 277,
%     pages = {1411-1429},
%    adsurl = {http://adsabs.harvard.edu/abs/1995MNRAS.277.1411U},
%   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
% }



\end{thebibliography}
\end{document}

